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An  analytical  solution  is  developed  to  describe  the  steady,  isothermal,  closed 
streamline  velocity  field  within  an  encapsulant  and  melt  in  a cylindrical  cavity  with 
uniformly  translating  walls  for  the  creeping  flow  limit  under  zero  gravity.  The 
solution  has  application  to  liquid  encapsulated  melt  zone  semi-conductor  growth. 

Equations  encompassing  both  no-slip  and  no-stress  boundary  conditions  are 
presented.  The  model  solutions  are  dependent  upon  three  parameters:  the  zone 
height,  the  unperturbed  interface  location  and  the  ratio  of  encapsulant  to  melt 
viscosity.  A linear  stability  analysis  is  applied  to  the  no-stress  solution  to  determine 
the  critical  value  of  surface  tension  at  which  the  melt  becomes  unstable.  Results 
obtained  by  varying  the  encapsulant  to  melt  viscosity  ratio  indicate  an  increase  in 
stability  as  this  ratio  is  lowered. 
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CHAPTER  1 
INTRODUCTION 


There  is  a significant  effort  currently  underway  to  grow  large  diameter 
Gallium-Arsenide  (GaAs)  semi-conductor  crystals.  One  method  being  considered  is 
that  of  a liquid  encapsulated  melt  zone  (LEMZ).  The  LEMZ  process  is  shown 
schematically  as  Figure  1.1.  In  this  process  a solid  rod  of  polycrystaline  material  is 
surrounded  by  a chemically  inert  encapsulant  and  the  entire  arrangement  is 
translated  relative  to  side  heaters.  The  polycrystaline  material  melts  and  re-solidifies 
as  a single  crystal.  The  LEMZ  method  is  being  examined  as  a way  of  minimizing 
diffusion  of  As  out  of  the  semi-conductor  melt,  thereby  resulting  in  a higher  quality 
crystal.  During  LEMZ,  convection  occurs  by  three  separate  processes:  1)  due  to 
density  gradients  (Rayleigh  effect),  2)  due  to  surface  tension  gradients  (Marangoni 
effect)  , and  3)  due  to  the  difference  in  the  solid  and  melt  densities  (Advective 
effect).  In  a zero-gravity  environment  the  Rayleigh  effect  can  be  ignored.  The 
Marangoni  effect  can  be  minimized  by  selecting  an  encapsulant  which  is  chemically 
inert,  but  of  much  higher  viscosity  that  the  melt.  Therefore,  it  is  anticipated  that  the 
largest  driving  force  for  convection  is  the  advective  effect.  In  order  to  study  the 
advective  effect,  surface  tension  will  be  considered  to  be  constant  throughout  this 
study,  thereby  eliminating  the  Marangoni  convection.  Since  only  the  advective  effect 
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FIGURE  1.1  Schematic  of  LEMZ  Process 
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is  considered,  a coordinate  transformation  to  a fixed  frame  of  reference  with  a 
moving  outer  wall  is  made.  This  model  system  is  shown  in  Figure  1.2.  The 
convection  pattern  for  this  system  is  closed  streamline  in  nature. 

Fluid  motions  with  closed  streamlines  generally  exhibit  two-  or  three- 
dimensional  flow  patterns.  The  case  of  cylindrical  geometry  with  axially  symmetric 
two-dimensional  flow  has  been  studied  for  a model  system  without  encapsulant.1  In 
this  case,  an  analytical  solution  was  developed  for  the  steady,  closed  streamline 
velocity  field  for  low  Reynolds  numbers.  In  this  creeping  flow  limit,  the  governing 
equations  are  linear  and  therefore  this  solution  is  unique.  The  method  of  solution 
used  is  applicable  to  the  LEMZ  model. 

The  stability  of  the  fluid  column  described  by  the  melt  is  a major 
consideration.  Large  diameter  crystals  of  long  length  are  desirable  from  the 
processing  viewpoint.  However,  it  has  long  been  known  that  a static  fluid  column 
becomes  unstable  and  breaks  apart  into  droplets  when  the  aspect  ratio,  length  by 
radius,  exceeds  the  circumference  of  the  column  (the  Rayleigh  limit).2-3  The  LEMZ 
model  has  two  complicating  factors,  the  closed  streamline  motion  of  the  melt  and  the 
encapsulating  fluid,  which  affect  this  criteria. 

The  stability  of  liquid  cylinders,  vertical  jets,  and  bridges  (liquids  between 
solid  surfaces)  has  been  studied  for  various  physical  properties  and  geometric 
configurations.  Many  of  these  systems  are  simpler  in  nature  than  the  model  under 
consideration  and  considerable  insight  can  be  gained  by  examining  them. 


4 


5 


The  simplest  case  is  that  of  an  inviscid  melt  and  inviscid  encapsulant  in  a 
infinite  coaxial  cylinder,  under  isothermal  conditions,  without  gravity,  subjected  to 
axi-symmetric  disturbances.  For  such  a system,  the  stability  is  independent  of  the 
surface  tension  and  adheres  to  the  Rayleigh  limit.4  Introducing  viscosity  into  the 
melt  phase  does  not  alter  this  result  and  further  indicates  that  the  stability  is  also 
independent  of  the  melt  viscosity.5  Recent  studies  of  this  system  have  subjected  it 
to  a "wind  stress"  shear  condition  at  the  melt-encapsulant  interface.6,7  It  was  found 
that  shear  can  lead  to  stability  beyond  the  Rayleigh  limit  for  a certain  "window"  of 
shear  values  far  beyond  the  creeping  flow  limit.  Therefore,  when  flow  is  present,  the 
stability  of  the  system  is  a function  of  the  melt  viscosity. 

Viscosity  may  also  be  introduced  in  the  encapsulant.  Doing  so  indicates  that 
the  viscosity  ratio  of  encapsulant  to  melt,  s,  determines  the  maximum  instability,  i.e. 
wavelength  of  the  fastest  growing  disturbance,  of  the  melt.8  Introducing  flow  results 
in  a class  of  problems  generally  referred  to  as  concentric  two-phase  pipe  flow.  Both 
the  effect  of  viscosity  and,  more  recently,  density  differences  on  stability  have  been 
examined.  It  has  been  shown  that  for  different  density  fluids  a stable  flow  can  always 
be  achieved  given  a sufficiently  large  pressure  driving  force.9  It  is  worth  mentioning 
that  if  the  flow  were  driven  by  a moving  wall,  rather  than  a pressure  drop,  that  the 
interface  could  not  be  vertical  in  the  presence  of  gravity. 

Closed  systems  have  not  been  studied  as  much  as  the  open  systems.  However, 
density  differences  for  an  inviscid,  stationary  closed  system  have  been  shown  to  affect 
the  growth  rate  of  disturbances.10  Studies  of  viscosity  effects  are  absent  in  the 
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literature.  The  ratio  of  encapsulant  viscosity  to  melt  viscosity,  s,  is  of  primary 
importance.  If  this  ratio  is  high  then  one  would  expect  the  encapsulant  to  begin  to 
act  like  a solid  and  increase  the  stability  of  the  melt.  However,  the  closed  streamline 
flow  which  occurs  in  the  encapsulant,  and  is  subsequently  transferred  to  the  melt, 
tends  to  de-stabilize  the  melt.  Thus  the  encapsulant  gives  rise  to  competing  effects, 
the  outcome  of  which  cannot  be  determined  a priori.  The  resolution  of  these  effects 
is  the  purpose  of  this  study. 


CHAPTER  2 

MODEL  FORMULATION 


The  dimensionless  equation  governing  the  flow  in  both  the  encapsulant  and 
melt  phases  is  (see  Appendix  A) 

L2T  = 0 (2.1) 

This  fourth  order  partial  differential  equation  can  be  decomposed  into  two  second 
order  partial  differential  equations  such  that 

LY  = rco  (2.2.1) 

L(a-q)  = 0 (2.2.2) 

7 is  obtained  by  first  solving  equation  (2.2.2)  for  o>  and  then  solving  equation  (2.2.1). 
The  boundary  conditions  for  equations  (2.2.1)  and  (2.2.2)  are  given  in  Appendix  C. 
The  no-flow  and  no-slip  velocity  conditions  given  in  Table  C.l  are  converted  to 
stream  function  conditions  by  using 
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Since  one  can  write  ijr'”  = T^-constant  = 0 and  Lij/71  = r o,  the  value  of  the  constant  can 
be  factored  out.  Therefore  it  can  arbitrarily  be  set  equal  to  zero.  Performing  similar 
calculations  for  the  other  10  boundary  conditions  in  Table  C.l  results  in  Table  2.1. 
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TABLE  2.1  Non-Interfacial  Y Boundary  Conditions 


MELT 

ENCAPSULANT 

rn  (r, 
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dxpm 
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z=0)  = 0 
= 0 

z =0 

re  (r, 

f y 
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=0,  z)  = 0 
= 0 

r=0 

1 

(r= 

/ 

aYe 

i ^ > 

1,  z)  = 0 
= 1 

r=l 

Since  the  boundary  conditions  for  each  phase  are  different,  each  phase  is  considered 
separately. 

Melt  Phase 

Solving  for  o 

The  boundary  conditions  for  o are,  in  general,  functions  of  z and  r,  i.e. 


g>(0,  z)  = 0 
a)(p,  z)  = kj?) 
a>(r,  0)  = kjr) 
<*>(r,  p)  = k3(r) 


The  value  of  <o(0,z)  is  known  due  to  the  velocity  boundary  conditions  of  the  melt  (see 
Table  C.l).  In  order  to  solve  equation  (2.2.2),  o>  is  written  as  a>  = CD1  + ci>2+t‘)3  where 
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the  boundary  conditions  of  Oj,  o>2  ar*d  ^>3  are  homogeneous  on  three  sides  and  are 
non-homogeneous  on  the  fourth  side  in  the  same  manner  as  g>.  This  results  in  three 
equations  of  the  form  L(r<y,)  = 0.  Each  of  these  equations  can  be  solved  by  the 
separation  of  variables  technique. 

Letting  o)J(r,z)=F/(r)G/(z)  and  substituting  into  equation  (2.2.2)  results  in 


1 d2G,- 

Gi  dz2 


1 

^7 


d2F,- 
r 

dr2 
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dF  / 

dr 


= Y / 


or 


+ r 


dF,- 

dr 


(2.3) 


and 


dG, 

—L  + G,Yf  = 0 (2.4) 

d z1 


Equation  (2.3)  may  be  reduced  to  the  familiar  "Bessels"  equation  by  using  an 
appropriate  variable  transformation.11  Thus  the  solution  to  equation  (2.3)  is 


if  y,->0 


(2.5.1) 


or 


F/  = ClIl(r/N~)  + C2Kl(r/ 1 Yj  I ) if  Y,  <0 


(2.5.2) 
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Equation  (2.4)  is  easily  solved  to  get12 


G,  = CjCoshjz^Y/--)  + C4sinh(z^y~) 


(2.6) 


The  constants  are  determined  by  applying  the  boundary  conditions. 

The  sign  of  Y/,  the  separation  constant,  is  determined  by  assuming  y/  to  be 
either  ( + ) or  (-)  and  determining  if  a contradiction  is  reached  when  the  boundary 
conditions  are  applied.  Using  as  an  example 


ui  = Fi  Gi 


subject  to 


o^jr,  0)  = Ojfr,  fi)  = QjjO,  z)  = 0 
cojfr?,  z)  = k^) 


1)  Assume  y>0 

2)  Apply  a)(r,  0)  = 0 

o = fJc^i)  * c4(0)] 

•C3  - 0 


3)  Apply  (o(r,fi)  = 0 


0 = F1C4sinh(/3v/y) 


.-.sinh(/3v^)  = 0 


12 


Pyfy  = 


nni 


nn ' 


P 


Which  contradicts  the  original  assumption  of  y>0.  Therefore,  y<0. 
4)  Apply  a >(r,  0)  = 0 


0 = 


FjC4sinh 


0 = F1C4sinh(/3//iy’j— ) 
-*.sinh(^iV  | Y I ) = 0 

1 3i\J  | y | = nni 


lYl  = 


nn 

T 


5)  Apply  o>1(0^)  = 0 


0 = [c,(0)  + CjMJGj 

■C,  = 0 


Thus 
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OO 

wi  = E where  an  = r^- 

n=  1 P 

Following  the  same  procedure  for  <i>2  and  o>3  results  in 

OO  r 

w2  = E JiM  + G^cosh^) 

/2  =1 
OO 


w3  = E C^Ji(v)sinhM 

/j  l 


where  bn  is  determined  from  J1(6/Jr1)=0.  Therefore, 


OO 


<0 


m 


n = 1 


= e c;;;  h(anr)  si4v) + fZl  sinh(M + cosh^)  ^(v)  (2.7) 


where  three  constants,  due  to  the  unknown  boundary  conditions,  remain  to  be 
evaluated.  Each  series  is  allowed  to  be  summed  from  1 to  °o  rather  than  0 to  ® 
because  a0=b0= 0 and  sin(fl^)=J1(6<?r)  = 0. 

Examination  of  equations  (2.3)  and  (2.4)  indicate  that  they  are  Sturm-Louiville 
in  nature.  Therefore,  the  characteristic  functions,  in  this  case  sin^^)  and  J i(bnr), 
which  arise  from  their  solution  are  each  a complete  basis  set.  This  allows  any 
function  to  be  expanded  in  an  infinite  series  of  these  characteristic  functions  and 
convergence  in  the  domain  is  assured.  Additionally,  due  to  the  Sturm-Louiville 
equations,  the  functions  are  also  orthogonal. 

Solving  for  Y 

Separating  7 into  two  parts  such  that 
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results  in 


Expanding  ¥ 


results  in 


L¥m  = L¥^  + L¥™ 
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E 0TlMsinM 

/i  1 


oo 


L¥?  = 


- E 

/i  =1 


C^sinhfv)  + C^cosh^jjrJ^) 


in  terms  of  sin(a/rz)  and  ¥™  in  terms  of  rJ1(6„r)  so  that 


*7  - E 


OO 
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OO  r 1 OO 

E sinM  = rE  Ij^rjsin^) 

n=l  /j=l 


(z)rJi(bnr)  = r£  sinhf^)  + cosh^)]  J^r) 

n=  1 


(2.8.1) 


(2.8.2) 


(2.9.1) 


(2.9.2) 
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Since  the  sin(a7rz)  are  independent  of  each  other,  as  are  the  rJ1(6/Ir),  every  n‘h  term 
of  equations  (2.9)  must  equal  each  other.  Taking  this  into  account,  and  performing 
the  indicated  operations,  results  in 


d7»  _ 1 dT>) 

dr2  r dr 


(2.10.1) 


d^(z) 

dz2 


bl  ^ sinh (*v)  + tZ  cosh (b„z) 


(2.10.2) 


Thus  equation  (2.2.1),  a second  order  partial  differential  equation  in  two  variables, 
has  been  converted  into  two  second  order  ordinary  differential  equations.  Equations 
(2.10)  can  be  solved  by  using  the  variation  of  parameters  technique.  The 
complimentary  solution  of  equation  (2.10.1)  is 


* Kt(v) 


and  a particular  solution  is 


( \ 

n = C*r 

_U 

ll(anr)  + rll(anr) 

D 

an 

'l  V 

thus 


O)  * CZrh(v)  * C^I^r) 


c 


m 
3 n 


(2.11) 
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The  complementary  solution  of  equation  (2.10.2)  is 


= C™  cosh^)  + CJ  sinh(v) 


and  a particular  solution  is 


= z 


C5n  coshM  + C"  sinh^JI  + -Ji  sinh(z>„z) 


thus 


T:£(z)  = (c£!  + cs”  cosh(^^)  + (c£  + C^zjsinh^)  (112) 


So,  finally, 


00  l 

= E[CSrIlM  + ^ ll(anr)  + CZ  r Kl(<v)]  sinM 


rt=l 

OO 

- E 

n=  1 


C*I  + C5KZ)COshM  + (Cto  + C7^Z)sinhM/Jl(6/Ir) 


(2.13) 


Encapsulant  Phase 

Solving  for  oi 

The  procedure  for  solving  for  0 in  the  encapsulant  is  the  same  as  for  the  melt. 
The  boundary  conditions  are 


Since  o in  non-homogeneous  on  all  boundaries,  o>  is  written  as  10  = q 3 + q-,  + g>3  + g>4 
and  four  constants  will  occur  in  the  final  equation  for  o.  Equations  (2.3)-(2.6)  are 
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o)(r?,  z)  = k-fe) 
o)(l,  z)  = k2(?) 
o(r,  0)  = k3(r) 

<*(r,  p)  = kjr) 

directly  applicable  in  the  encapsulant.  Their  solution  results  in 


OO 


0)^ 


- 

n =1 


tanh(c^)cosh(c;rz)  - sinh(c^)  ] P^r) 


“4 = e ci,sinhMpiM 


«=i 


thus 


(2.14) 


where  cn  is  determined  from  P1(c/I)  = 0 and 
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Once  again,  since  c0  = 0 and  P1(c^r)  = 0,  each  series  is  summed  from  1 to 
Solving  for  7 

The  characteristic  functions  from  the  eigenvalue  expansion  in  the  encapsulant 
are  sin(a/rz)  and  P^c^r).  As  was  the  case  for  the  melt  phase,  these  functions  are 
complete  and  orthogonal  (See  Appendix  D for  the  proof  of  the  orthogonality  of  the 
P^s).  The  solution  procedure  for  T¥e  is  analogous  to  that  for  r¥m.  Thus 


EL 

/I  =1 


0° 

Yl*(r)  si4v)  = E r Ce\n  KiM  + C^Ii (anr)  sin (altz)  (2-15) 


n=  1 


oo 


E 

n=  1 


(2.16) 


Equating  terms  of  the  orthogonal  functions  and  solving  the  resultant  ordinary 
differential  equations,  again  by  using  variation  of  parameters,  results  in 


Ct,  IiM 
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and 


^l/|(r)  ^1/1  r h(anr)  + ^2 n r ^1  [anr)  + ^3n  ^o[anr)  + ^4n  r 2 ^o[c 


v)  <2-17) 


(2.18) 


CHAPTER  3 

LINEAR  STABILITY  ANALYSIS 


A linear  stability  analysis  is  performed  to  determine  the  point  at  which  the 
interface  is  no  longer  vertical.  A deflection  from  verticality  is  considered  unstable. 
The  stability  of  the  interface  is  determined  by  expanding  rj  and  functions  of  p in 
MacClaurin  series  about  an  infinitesimal  perturbation,  e,  and  then  linearizing  the 
series  by  dropping  terms  of  0(e2)  and  higher*.  In  general 

/(e)  = m + ef'(0)  ♦ - * e'‘  ^('‘^(0)  * K 
thus,  in  linear  form, 


The  perturbation  parameter  e is  not  defined  in  this  study,  as  e is  not  necessary 
in  the  linear  analysis  in  order  to  determine  the  value  of  Ca”1.  Should  one  wish  to 
determine  the  branching  into  the  new  steady  solution,  the  non-linear  problem  of 
0(e)  must  be  solved,  at  which  time  e must  be  defined.  That  we  can  determine 
exactly  the  point  at  which  this  branching  occurs  from  the  linear  problem  requires 
some  consideration.  We  shall  see  that  the  linear  stability  analysis  results  in  the 
eigenvalue  problem  Cp  = Cac1p.  It  should  be  realized  that  when  the  original  non- 
linear equations  are  perturbed  and  linearized  that  the  eigenvalue  must  also  be 
perturbed  and  linearized  so  that  Ca'1(e)  = Ca^1(0)  + eCa-1'(0).  It  is  the  O(e0) 
term  of  the  eigenvalue  which  appears  in  the  final  eigenvalue  problem.  It  is  apparent 
that  this  term  is  independent  of  e and  that  the  linearization  does  not  affect  its  value. 
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r?(e)  = p(0)  + €T\\  0) 

/(b,<0  = /(b,e)|t=0  + e 


d/(r?,e) 


de 


(3.1) 


e=0 


where/ may  be  a variable,  a function,  or  a functional.  Since  p is  itself  a function  of 
e,  the  chain  rule  must  be  applied  to /such  that 


/(b,e)  = /(p,e)  |e=0  + e 


= /(h,e)  |e=0  + e 


de  + ffip.e)  dr? 
, 3e  de  6lr  drj  de 

' df{r],e)  dp  df(p,e)  dr 


le=0 


de 


de  6^  dr] 


e=0 


= f + * 

= f + £/[0] 


f ~cf  dr 

/(0)  + ”(0)1^ 


This  final  form  introduces  new  notation  as  indicated.  Writing  p in  this  same  notation 
results  in 

n = ri  + ep(0) 


It  should  be  noted  that  except  along  the  interface,  where 


!-/(0)-/[0]- 


Perturbing  L2T=0  in  this  manner  results  in  two  separate  problems,  L2T  =0, 
the  base  problem,  and  L27j0pO,  the  perturbed  problem  (see  Appendix  E for  a 
detailed  discussion  of  the  perturbation  technique  resulting  in  these  equations).  To 
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solve  these  equations  the  boundary  conditions  must  also  be  perturbed.  For  the  base 
problem,  the  conditions  listed  in  Table  2.1  hold  with  ¥ replaced  with  7.  For  the 
perturbed  problem  the  conditions  listed  in  Table  2.1  hold  with  ¥ replaced  with  7(0), 
except  for  the  Neumann  condition  at  r=  1 which  must  equal  0.  The  interfacial 
conditions  are  considered  separately. 

The  solution  procedure  for  these  equations  is: 


1)  Solve  L27  = 0 for  ¥m  and  ¥e. 

2)  Solve  L“T[o]  for  ^O)  anc*  ^(0)  *n  terms  of  b(o)  and  7. 

3)  Substitute  7^  into  the  normal  component  of  the  momentum 
balance  and  solve  the  resulting  eigenvalue  problem  for  rj,0y 


Kinematic  Conditions 


Base  Problem 


In  the  base  problem,  rj=r j.  Therefore  equation  (C.2)  reduces  to 


(3.3.1) 


and  equation  (C.3)  reduces  to 


(3.3.2) 


Perturbed  Problem 


Perturbing  equation  (C.2)  results  in 


23 


- e 


f 

\ 

/ 

dr]  dYm 

dYm 

= 4-  £ 

dYm 

dr  J 

[0]  & 

{ dz  J 

[0] 


equating  0(e)  terms  and  expanding  the  terms  yields 


Y 


m 

(0) 


m 


r,  ’ ' ”«»  "IT 


» 


m 

(0) 


(3.4.1) 


Using  the  same  procedure  for  the  encapsulant  results  in 


Y 


(0) 


"(0) 


dY‘ 


(0) 


(3.4.2) 


No  Slip 


Base  Problem 

Equation  (C.6)  reduces  to 


/ > 

( \ 

3Y/71 

CL> 

a 

i dr  J 

r. 

l dr  J 

(3.5.1) 


Perturbed  Problem 

Using  equation  (C.6)  and  following  the  same  procedure  as  above  results  in, 


/ ™ \ 

dY 

d*(0) 


m 


g(o  y 


/ e \ 

!!w 

dr  J 


+ *(0) 


(3.5.2) 
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where 


m _ 

S(o)  - h(0) 


a2?"1 

T2"" 


. e 

and  *(0)  * 1(0)  — - 
or 


Base  Problem 


Tangential  Stress 


Equation  (C.10)  reduces  to 


dZxPm 

~d^~ 


1 dWm 
+ 

rj_  dr 


\ 


J r. 


= S 


£xpe 

~d^ 


\ 


i a?e 

r±  dr 


(3.6.1) 


Perturbed  Problem 


Using  equation  (C.ll)  and  following  the  above  procedure  results  in 


V> + 


a2? 


<2>  + 


a? 


rj  dr 


gxvm 

(0)  + <^*(0) 


dzx 


(3.6.2) 


l(0)  + 


*%)  . 1 ay?0)  , 
dr 2 *1 


a2? 


(0) 


dzA 


where 


i m _ 

>)  - 


4a» 


dz 


/ _ 

a2?"1 

drdz 


17  (0) 


a3Ym  a3Tm  i aY™ 


cKfe"1 


ar- 


2 ar 


-'/■i 


and 
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= 4 


dr\ 


(0) 


dz 


/ _ \ 

a2ye 

cfo2r  > 


+ n 


(0) 


a3?* 

c^ck'1 


a3Te 

&-3 


i ay* 

_2  dr 


Normal  Component  of  Momentum 


Base  Problem 

Equation  (C.14)  reduces  to 


(pm 


- Pe 


Ca"1  2 

f 

a2?'” 

> 

gTpe 

lr i - = - 

1 ri  h 

drdz 

- s 

drdz  , 

(3.7.1) 


This  equation  is  not  required  for  the  evaluation  of  the  constants.  However,  it  would 
be  required  to  determine  pressure  differences. 

Perturbed  Problem 

Using  equations  (C.14)  and  (3.7.1)  and  following  the  above  procedure  results 
in 


26 


[0]  1 [0]/ri 


m-PL ) +Ca_1 


'(0) 


a2 


c fz“ 


i aV” 

1 ^*(0) 

r1  arciz 


1 dxp(0) 
2 dz 


- s 


1 


^*(0) 


1 ^(0) 


rl 

drdz 

2 

rl 

dz 

11  (0) 

a3^"1 

3 

a2*?"1 

rl 

dr2dz 

ri 

drdz 

^(0) 

3 

rl 

dr2dz 

ri 

drdz 

(3.7.2) 


An  examination  of  equation  (3.7.2)  reveals  that  7/?I  and  are  known  from  the  base 
problem  solution  and  q)  ar*d  ^(o)  are  known  in  terms  of  from  the  perturbed 
problem  solution.  Thus,  the  right-hand-side  of  this  equation  is  solely  a function  of 
P(0).  In  order  to  obtain  the  desired  eigenvalue  problem,  the  perturbed  pressure 
difference  must  also  be  known  in  terms  of  rj^y  This  difference  is  derived  in 
Appendix  E.  Expanding  rj^  in  terms  of  the  appropriate  sine  functions  results  in  the 
eigenvalue  problem 

( A - Ca-1Z?)x  = 0 (3,8) 


and  the  non-trivial  solution  of  equation  (3.8)  is  determined  by 


= 0 


DET 


A - Ca ~lB 


(3.9) 
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As  we  shall  see  later,  for  aspect  ratios  beneath  the  Rayleigh  limit  the  largest  real 
eigenvalue,  Ca'1,  determines  the  point  at  which  the  system  becomes  unstable.  Since 
the  characteristic  viscosity  and  velocity  are  fixed,  this  in  turn  determines  the  critical 
value  of  surface  tension  for  the  system  to  be  stable.  If  only  complex  eigenvalues 
result  from  the  solution  of  equation  (3.8)  the  system  is  oscillatory,  which  is 
considered  unstable  in  this  study. 


CHAPTER  4 
MODEL  SOLUTION 


Equation  (2.13)  defines  the  stream  function  for  the  base  problem  melt  phase. 
The  required  seven  boundary  conditions  are  given  by  equations  (3.3.1),  (3.5.1)  and 
Table  2.1  (Although  this  adds  up  to  eight  boundary  conditions,  it  should  be  noted 
that  the  condition  for  bounded  centerline  velocity  has  been  used  in  determining  <i>.). 
Applying  the  Dirichlet  conditions  results  in 


(4.1) 


where 
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yZ  = sinM 

yZ  = r]\[bnr ) 


Equation  (2.18)  defines  the  stream  function  for  the  base  problem  encapsulant 
phase.  The  required  eight  boundary  conditions  are  given  by  equations  (3.3.2),  (3.6.1) 
and  Table  2.1.  Applying  the  Dirichlet  conditions  results  in 


OO 


- E 

/!=1 


€ y~i€  \ € 

c \F- 


2n]y\n 


r,e 

C3/A3/i 


€ j-i€  \ € 

4/i  4n)y2n 


(4.2) 


where 


r\ M + B^rh 

M 

1 + B3nrKl[anr) 

H 

!«-.) 

r2Manr) 

1 + B^r\A 

.w 

1 ♦ B5nr K,| 

M 

Vl) 

Fl  = j(P~z)  sinhfc^) 


-2zsinh 


■Jp 
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and 


B 


2 n 


. rManrl)Kl(an)  ~ Kl{anrl)lo(an) 


B 


bt 


R _ ll(anrl)lo(an)  ~ rlll(an)lo(anrl) 

°3n  n 

Bbi 

R _ rlKb(anrl)Kl(an)  ~ 

4 n B 

°l/i 

R _ Ii(e*«ri)K0(a«)  - rManri)h{an) 

5n  ' 

= Kl(«n,'l)Il(fl«)  - IiK''i)KiK) 

It  is  easily  seen  that  equations  (4.1)  and  (4.2)  satisfy  the  Dirichlet  conditions.  The 
Neumann  conditions  are  applied  by  expanding  the  z-functions  in  fourier  sine  series 
and  the  r-functions  in  fourier  bessel  series  (rJ1(h/1r)  for  the  melt  phase  and  rV^cj) 
for  the  encapsulant  phase).  The  boundary  conditions  are  applied  in  the  following 


order: 
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1)  melt  phase  at  z = 0 to  eliminate  C^n  in  terms  of  C?%n  and  Cf^n. 

2)  melt  phase  at  z=/3  to  eliminate  C%n  in  terms  of  C^n. 

3)  encapsulant  phase  atz  = 0 to  eliminate  in  terms  of  C^,,  and  C^,. 

4)  encapsulant  phase  at  z=(3  to  eliminate  C^,  in  terms  of  and  C?hl. 

5)  no  slip  at  the  interface  to  get  C^n  in  terms  of  C^,,  and  C^,. 

6)  continuity  of  tangential  stress  at  the  interface  to  get  in  terms  of  C^,, 
and  C^. 

7)  use  the  resultant  equations  from  (5)  and  (6)  to  eliminate  and 
obtain  in  terms  of  C^. 

8)  encapsulant  phase  at  r=l  to  get  C^,  in  terms  of  the  fourier  sine 
coefficients  of  1. 

This  results  in 

_ OO  OO  OO  D TIf  r 

ne  V'  V nk^kr  lm  ‘•m  - ... 

cs,  = -L  E E — (4-3.1) 


where  the  are  the  fourier  sine  coefficients  of  1.  Rnk  is  the  inverse  of  a matrix 
which  arises  from  step  (8)  above.  It  is  defined  by 


OO 


i =i 


+ E cjMik 

i =i 


= 8 


nk 


W/m  the  inverse  of  a matrix  which  arises  from  step  (7)  above.  It  is  defined  by 
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where 


OO 

E wik 

k=  1 


O 


5/t/i  + — 


btk 


^'(1) 


= <S 


/« 


R = n Flk  (rl)  ^ n °4kp 

Bik  - QVk  ^7T7  ' lkSjk  + E 01;>  TrT 


-1 


C = O (ri 

^ “ CV'  ^7 
Fi i h 


a. 


- T2lsjl  + E Sl;>  — 319 


p=L 


Mlk  - 


e ; _1  J r^e 


E 

fJt'O)  1-1  ^'(D 


*s  inverse  of  matrix  which  arises  from  step  (5)  above  and  is  defined  by 


E Qyk 

k=l 


Jkn 


°5nk 


Tlk  ~ 


li{akri)Kiakri) 
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OO 


E iL(V“nfUfLPCn  +flp>- 


m =1 


e 

mn 


oo 

^2 np  = ^2  ^'bn^>anf'2nm{7>mp^mn  + f 4mp  *■ 
m =1 


e 

mn 


OO  _ r_ 

®3np  = Y1  ^7m  [r l)  a n f \nm\f  3mn  ^ mn  + f 4mn^> 


m =1 


3 mp  ^ mn  J 4mp  mn 


oo 


°mp  ’ E yLh WflJfLpiL  +fLp\ 


m=  1 


> 3mp  ^ mn  J 4mp  mn 


(-)-  = V'  Vm  I ) f If"1  rm  f"!  ; "J 

5np  — ^ ^2m  \ 1/  anJ  Inm^  3mp  ^mn  + J 4/np  mn 
m =1 

and  5 is  the  kronecker  delta,  ~j\nm  are  the  fourier  P1  coefficients  of  F\ n,f2/im  are  the 
fourier  Px  coefficients  of  F^,  ~jylm  are  the  fourier  sine  coefficients  of  F%n^Anm  are 
the  fourier  sine  coefficients  of  P%t,finm  are  the  fourier  Jx  coefficients  of  F7„»/3/im 
are  the  fourier  sine  coefficients  of  F"ln,  ~f\ntn  are  the  fourier  sine  coefficients  of  F^n 
(see  Appendix  F for  all  fourier  coefficients)  and 

tc  . i Mr1 

cmP  ~ (-lfsinh (cmJ3) 
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re  _ _ 1 P 

^ mn  ; 7 

cm/3  - (-l)nsinh[cmp) 

~m  _ 1 (-I)"*1 

Kmn  ^ ; r 

2 bmp  - (-lfsinh (bmp) 


1  P 

2 bmP  ~ (-Ifsinh (bmp) 


These  O matrices  arise  from  the  elimination  of  the  constants  associated  with  the  z- 
functions  in  steps  (1-4)  above. 


OO 


= E Wkn 

k = 1 


OO 

E C2m°2/nk 

m=  1 


(4.3.2) 


OO 


e 

nm 


oo  _ r_ 

E am  ^ nm  ^1 rJ  lmn 

m =1 


= r 


Ce 

ln^ln 


+ 


T2n 


(4.3.3) 


(4.3.4) 


(4.3.5) 
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OO 

CZ  - E 

m=  1 


(4.3.6) 


OO 

p/ti  ^ i "i  p»n 

'“4/1  ~~  ^ am^nm  1 nv  \mn 
m=  1 


(4.3.7) 


A close  examination  of  these  equations  reveals  that  if  n is  an  even  number  then 

q*=c^=c7„=o- 

Equation  (2.13)  also  defines  the  stream  function  for  the  perturbed  problem 
melt  phase.  The  required  seven  boundary  conditions  are  given  by  equations  (3.4.1), 
(3.5.2)  and  Table  2.1.  Applying  the  Dirichlet  conditions  results  in 


¥ 


m 


n=  1 


(Ofsn  ~ CbtFbt)^bi 


^ m r m 
Sr  3n 


^ m ~ m \-  m 
C4nF4n  K, 


(4.4) 


where 


Fm  - P T7m 

^3 n - 7 rF3/i 

2sinh(6,J£j 


-1 

2cosh|6/(/3j 
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and  fn  (0)  are  fourier  sine  coefficients  of 


dx¥m 

"(0)— 


Equation  (2.18)  defines  the  stream  function  for  the  perturbed  problem 
encapsulant  phase.  The  required  eight  boundary  conditions  are  given  by  equations 
(3.4.2),  (3.6.2),  and  Table  2.1.  Applying  the  Dirichlet  conditions  results  in 


oo 


'£  fn 
n =1 


(0)  5n 


Ce  Fe 

^lnr\n 


e 


4 nr2n 


(4.5) 


where 


P Fe 

7 \ 3« 

2sinh  (c„j9) 


-1 

2cosh|c/I/?J 


h 

K)Ki 

V)  ’ h 

v)K  ia> 

.)! 

l'l( 

an)Kl(anr l)  - ri 

an\ 

K i(a„) 

and^(0)  are  the  fourier  sine  coefficients  of  n . ^ & . 

(0)  dr 

A strong  analogy  to  the  base  problem  is  seen  in  these  equations.  The 
additional  functions  of  r arise  from  the  kinematic  conditions.  All  other  functions 
could  be  equated  to  the  base  problem  functions  through  an  appropriate  adjustment 


in  the  constants. 
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It  is  easily  verified  that  equations  (4.4)  and  (4.5)  satisfy  the  Dirichlet 
conditions.  The  Neumann  conditions  are  applied  in  the  same  manner  as  in  the  base 
problem.  The  resultant  equations  for  the  constants,  in  terms  of  coefficients  of  the 
eigenvector  expansion,  are 

OO 

c£  = E°>,(o, 

1=1 


OO 

c£=I>S>/(0) 

/ =1 


c;  = eo,(o)  <4-6-3> 

1=1 


oo 

c£  - E (4-6-4) 


OO 

cl-E^TO  <4-6-5> 


Xj  a3 til71 1(0 ) 


(4.6.6) 


OO 

£ a4 nfl  1(0 ) 


(4.6.7) 


where 
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a2 nl 


OO  oo 

E Rnk?kl  + *nk  "Ifsk  (!)  + E^iO) 
*=1  r=l 


alnl 


OO 


-E- 


/i/t 


‘-1 


OO 


“3/1/ 


OO 


E* 


A:  =1 


e 

•n* 


w 


e r e 
1 kf  5kn 


e r£  e j-e 

a \kU  \kn  + a2klJ 2kn 


a4 nl 


OO 


E»e  e r e 

; nk  W\kV  5kn 
k= 1 


e fe 
a\klJ\kn 


e re 
a2klf  2kn 


m 

al nl 


w 


m 

lnl 


m e 

+ v2 nl  ~ ™2 nl 


T1 na\nl 


+ T2na2nl 


oo 


*=1 


m rfn 
a\klJ\kn 


oo 


k=  1 


and 
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OO  OO 


OO  ft  r U7 


Hi  - E £ 


n=l  m=l 


In  these  equations,  the  w;’s  arise  from  the  kinematic  conditions  and  are  the 
coefficients  of  the  eigenvector  expansion  of  the  fourier  sine  coefficients  of  th 
the  v^’s  arise  from  the  continuity  of  tangential  stress  at  the  interface  and  are  the 
coefficients  of  the  eigenvector  expansion  of  the  fourier  sine  coefficients  of  the  /z^™. 
The  ekl  term  arises  from  the  no  slip  condition  at  the  interface  and  is  given  by 


The  w2’s  are  the  coefficients  of  the  eigenvector  expansion  of  the  fourier  sine 
coefficients  of  the  g^oj- 

These  coefficients  are  given  by 


e m 
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Three  additional  "O"  matrices  arise  because  of  the  kinematic  conditions.  They  are 
defined  by 
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At  this  point  the  eigenmatrix  can  be  formulated.  Substituting  equations 
(E.12),  (4.1),  (4.2),  (4.4),  and  (4.5)  into  equation  (3.7.2)  defines  .4  and  B.  They  are 
given  by 
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which  arises  from  the  perturbed  melt  pressure  term, 
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which  arises  from  the  perturbed  stream  function  terms, 


which  arises  from  the  base  stream  function  terms.  The  encapsulant  is  analogous  to 
the  melt.  Thus, 
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CHAPTER  5 

DISCUSSION  AND  RESULTS 


The  base  problem  of  a vertical  interface  is  a solution  for  this  system.  In  the 
perturbed  problem  this  interface  has  become  deflected.  The  driving  force  for  this 
deflection  is  the  normal  stress  at  the  interface.  This  is  easily  seen  by  examining 
equation  (3.7.1).  Rewriting  this  equation  as 

AP  = Ca_1{x}  + {y}  (5-1) 

reveals  that  the  pressure  difference  between  the  melt  and  encapsulant  is  balanced  by 
the  surface  tension  and  the  normal  stress.  For  a fixed  value  of  surface  tension,  the 
pressure  difference  increases  as  the  normal  stress  increases.  At  some  point,  an 
increase  in  the  normal  stress  can  no  longer  be  sustained  by  the  pressure  difference. 
When  this  occurs,  the  interface  will  deform  from  the  vertical  to  achieve  a 
configuration  with  lower  surface  energy,  thereby  allowing  the  pressure  difference  to 
once  again  sustain  an  increase  in  the  normal  stress.  For  incompressible  fluids  the 
form  of  this  deflection  must  be  such  that  the  total  volume  is  conserved.  Therefore, 
only  deflections  as  shown  by  figure  5.1-c  are  possible.  This,  of  course,  assumes  that 
the  liquids  continue  to  completely  wet  the  bounding  solids.  It  also  discounts  the 
possibility  of  the  solid-melt-encapsulant  line  shifting.  The  deflections  depicted  by 
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figures  5.1-a,b  are  only  possible  if  the  distance  between  the  upper  and  lower  plates 
changes,  which  may  occur  in  an  experiment  where  there  is  a vapor  space  above  the 
encapsulant.  The  point  at  which  the  deformation  will  occur  cannot  be  predicted 
quantitatively  from  equation  (5.1). 

Once  Ca”1  is  determined,  one  must  be  able  to  say  whether  it  is  to  the  left  or 
to  the  right  of  Ca”1  that  stability  is  found.  To  determine  this  we  examine  the  no 
flow  problem  and  then  extend  it  to  the  flow  problem  using  continuity  of  solutions. 
For  a viscous  liquid  jet,  Chandrasekhar5  finds  that 

q <*  -Ca(,1x  (5-2) 


where 


x = 1 - 


t 2 

hr. 


P 


(5.3) 


and  q is  a growth  constant.  For  aspect  ratios  below  the  Rayleigh  limit,  x is  negative, 
thus  q is  negative  and  disturbances  decay.  Suppose,  however,  that  one  fixes  q to  be 
some  negative  number  and  determines  Ca~ 1 then  makes  q more  positive  and  again 
determines  Cac  1 and  so  on  until  q is  positive.  At  this  point  Ca~ 1 becomes  negative, 
a physically  meaningless  situation,  but  one  which  serves  to  point  out  that  for  the 
viscous  liquid  jet,  decreasing  Ca“ 1 approaches  the  region  of  instability.  This  indicates 
that  the  region  of  stability  is  to  the  right  of  Ca  c 1 for  aspect  ratios  below  the  Rayleigh 
limit.  Similar  reasoning  leads  to  the  conclusion  that  aspect  ratios  above  the  Rayleigh 
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FIGURE  5.1  Typical  Deformations  of  Cylinders,  a)  Deformation  resulting  from 
compressing  the  zone;  b)  Deformation  resulting  from  stretching  the 
zone;  c)  Possible  deformation  of  the  interface. 
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limit  result  in  increasing  Ca”1  approaching  the  region  of  instability  and  the  region 
of  stability  is  to  the  left  of  Ca”1.  This  behavior  should  also  hold  true  in  the  model 
under  consideration.  This  is  because  a small  flow  should  not  result  in  a dramatic 
reversal  of  the  stable  and  unstable  regions  of  the  no  flow  case. 

The  quantitative  determination  of  ac  is  made  by  using  equation  (3.9).  The 
order  of  the  eigenmatrix  must  be  sufficient  to  result  in  a converged  eigenvalue,  as 
that  is  the  quantity  of  interest.  However,  it  is  unlikely  that  this  order  will  result  in 
an  eigenvector  which  satisfies  the  normal  component  of  the  momentum  balance13 
(equation  3.7.2).  Additionally,  the  order  of  the  eigenmatrix  must  be  even.  This  is 
apparent  if  one  considers  the  movement  of  the  outer  wall.  Suppose  wall  movement 
downward  in  the  ( + ) z direction  with  velocity  + vz  results  in  Ca~3  for  a given  rl  and 
P-  Now  consider  wall  movement  upward  with  velocity  -vz.  The  only  change  will  be 
that  Cac1  will  also  change  signs.  Thus  Ca'1  will  always  occur  in  ± pairs.  An 
eigenmatrix  of  odd  order  must  possess  at  least  one  unpaired  eigenvalue  (Unless  the 
unpaired  eigenvalue  equals  zero,  which  is  a ± pair  unto  itself.  However,  a zero 
eigenvalue  indicates  a singular  eigenmatrix  which  cannot  occur  from  the  linear 
equations.)  and  therefore  is  disallowed. 

Although  the  order  of  the  eigenmatrix  must  be  sufficient  to  converge  the 
eigenvalue,  what  is  even  more  critical  are  its  elements,  each  of  which  is  an  infinite 
sum  unto  itself.  Thus,  for  a given  order,  each  element  of  the  matrix  must  also  be 
converged.  In  fact,  such  summations  appear  throughout  the  problem  (for  example 
the  O matrices)  and  in  each  case  convergence  must  occur  prior  to  using  the  quantity 
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in  the  next  step.  This  results  in  non-square  matrices  everywhere  except  for  the 
eigenmatrix  itself.  For  example,  w^kl  is  non-square  with  an  /-dimension  equal  to  the 
order  of  the  eigenmatrix  and  a ^-dimension  sufficient  to  result  in  a converged  a ^nl. 
The  extent  of  the  summations  required  for  convergence  varied  between  different 
quantities.  In  each  case,  convergence  was  assumed  when  the  value  of  the  summation 
remained  unchanged  upon  increasing  the  number  of  terms  in  the  summation  by  at 
least  50%  (for  example  from  200  to  300  terms).  Convergence  of  these  various 
summations,  particularly  in  the  perturbed  problem,  is  a non-trivial  matter. 

In  the  base  problem,  it  was  found  that  the  stream  functions  in  the  encapsulant 
and  melt  converged  for  50  terms  in  the  final  summation  (equations  (4.1)  and  (4.2)). 
At  the  same  time,  the  inner  summations,  in  particular  the  O matrices,  required  100 
terms.  These  matrices  require  additional  terms  because  they  result  from  the  fourier 
expansions.  These  fourier  expansions  expand  exponential  functions  in  terms  of 
periodic  functions  resulting  in  slow  convergence.  The  final  proof  of  convergence  is 
the  satisfaction  of  the  boundary  conditions.  Since  the  base  problem  is  linear  in 
nature,  there  can  be  only  one  set  of  constants  which  satisfy  the  Neumann  conditions 
(the  Dirichlet  conditions  are  automatically  satisfied  from  the  form  of  the  stream 
function).  Tables  5. 1-5.5  show  the  satisfaction  of  the  boundary  conditions  for  typical 
parameters.  The  degree  of  satisfaction  of  these  conditions  does  vary  somewhat  with 
the  value  of  the  parameters.  However,  it  was  found  that  the  50/100  term 
summations  achieved  an  acceptable  level  of  satisfaction  for  all  parameter  sets 
investigated. 
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In  the  perturbed  problem,  additional  series  arise  due  to  the  dependence  on 
the  base  problem.  Thus,  three  separate  levels  of  summation  are  required,  1)  due 
to  the  order  of  the  eigenmatrix,  2)  due  to  the  base  problem  dependence,  and  3)  due 
to  fourier  expansions. 

The  Ofyr^)  and  matrices  were  found  to  require  much  higher  levels  of 

summation  than  the  O matrices  from  the  base  problem.  Whereas  in  the  base 
problem  100  terms  were  sufficient,  300  terms  were  required  in  the  perturbed 
problem.  This  is  directly  associated  with  the  fourier  expansion  of  the  F5  functions, 
which  were  found  to  be  very  slowly  convergent. 
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TABLE  5.1  Melt  no-slip  boundary  condition  for  ^ = .5,  0 = .5,  5 = 3.  . 


r 

dxFm 

e\ 

raluated  at 

dz 

z = 0 

z=fi 

0.000 

0. 

0. 

0.025 

1.536  TO'9 

1.152  TO*4 

0.050 

1.369  TO'9 

6.060  TO'4 

0.075 

-3.053  TO'9 

13.607  TO'4 

0.100 

-1.223  TO'9 

24.111  TO'4 

0.125 

4.454  TO'9 

37.489  TO'4 

0.150 

0.53  MO'9 

53.603  TO'4 

0.175 

-5.816-10'9 

72.236  TO'4 

0.200 

0.816T0'9 

93.056  TO'4 

0.225 

7.078  TO'9 

115.571  TO'4 

0.250 

-3. 152  TO'9 

139.073  TO'4 

0.275 

-8.544  TO'9 

162.567  TO'4 

0.300 

6. 182  TO'9 

184.692  TO'4 

0.325 

10.276  TO'9 

203.618  TO-4 

0.350 

-11.004  TO'9 

216.954  TO'4 

0.375 

-12.063  TO*9 

22 1.682  TO'4 

0.400 

19.352  TO'9 

214. 160  TO'4 

0.425 

-5.309  TO'9 

190.340  TO'4 

0.450 

-269.877  TO'9 

146.493  TO'4 

0.475 

-2097. 149  TO'9 

81.160T0'4 

0.500 

0. 

0. 
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TABLE  5.2  Encapsulant  no-slip  boundary  condition  for  ^ = .5,  /3  = .5,  5 = 3.  . 


r 

evaluated  at 
dz 

z = 0 

0.500 

0. 

0. 

0.525 

-10.849T0'4 

-10.849  TO'4 

0.550 

1.793  TO'4 

6.875  TO'4 

0.575 

11.333-10"4 

8.067  TO'4 

0.600 

-3.845  TO'4 

-13. 131  TO'4 

0.625 

-11.788  TO'4 

-1.872  TO'4 

0.650 

6.275  TO'4 

17.374  TO'4 

0.675 

12. 194  TO'4 

7.418  TO'4 

0.700 

-9.286  TO'4 

-18. 144  TO'4 

0.725 

-12.540T0'4 

19. 142  TO'4 

0.750 

13.227  TO'4 

13.911  TO'4 

0.775 

12.780  TO'4 

-32.349  TO'4 

0.800 

-18.794  TO'4 

-2.822  TO'4 

0.825 

-12.801  TO'4 

45.984  TO'4 

0.850 

27.598  TO'4 

-18.268  TO'4 

0.875 

12.279  TO'4 

-59.345  TO'4 

0.900 

-43.786  TO'4 

57.217T0"4 

0.925 

-3.454  TO'4 

68.205  TO'4 

0.950 

143. 199  TO'4 

-207.086  TO'4 

0.975 

432.537  TO'4 

-591.893  TO'4 

1.000 

0. 

0. 
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TABLE  5.3  Moving  wall  no-slip  boundary  condition  for  ^ = .5,  /3  = .5,  s = 3.  . 


z 

dV e 
* 1=1 

0.000 

0. 

0.025 

0.990 

0.050 

1.048 

0.075 

1.001 

0.100 

0.975 

0.125 

1.000 

0.150 

1.018 

0.175 

1.000 

0.200 

0.984 

0.225 

1.000 

0.250 

1.015 

0.275 

1.000 

0.300 

0.984 

0.325 

1.000 

0.350 

1.018 

0.375 

1.000 

0.400 

0.975 

0.425 

1.001 

0.450 

1.048 

0.475 

0.988 

0.500 

0. 
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TABLE  5.4  Interfacial  no-slip  boundary  condition  for  ^ = .5,  /3  = .5,  s= 3.  . 


z 

dxpm 

I » L, 

dVe 

i » L, 

0.000 

0. 

0. 

0.025 

-1.133  *10'3 

-1.298 -10’3 

0.050 

-4.339 -10'3 

-4.340 -10'3 

0.075 

-9.183-10'3 

-9. 183  TO'3 

0.100 

-15.112-10'3 

-15.114-10'3 

0.125 

-21.506 -10'3 

-21.506  TO'3 

0.150 

-27.744 -10‘3 

-27.743  TO'3 

0.175 

-33.254 -10'3 

-33.254  TO'3 

0.200 

-37.556 -10’3 

-37.556  TO'3 

0.225 

-40.288 -10'3 

-40.288  TO'3 

0.250 

-41.225 -10'3 

-41.224  TO'3 

0.275 

-40.288 -10'3 

-40.288  TO'3 

0.300 

-37.556 -10*3 

-37.557T0'3 

0.325 

-33.254 -10'3 

-33.254  TO'3 

0.350 

-27.744 -10'3 

-27.743  TO'3 

0.375 

-21.506 -10‘3 

-21.506  TO'3 

0.400 

-15.112-10'3 

-15. 114  TO'3 

0.425 

-9.183-10"3 

-9.183  TO'3 

0.450 

-4.339 -10'3 

-4.341  TO'3 

0.475 

-1.133  -10‘3 

-1.187-10'3 

0.500 

0. 

0. 
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TABLE  5.5  Interfacial  tangential  stress  boundary  condition  for  ^ = .5,  /3  = .5,  s= 3.  . 


2 

d2xPm  x 1 dWm 
dr 2 r & 

J'='l 

d 2Ve  ldWe 

s + 

dr 2 r & 

0.000 

0. 

0. 

0.025 

-0.014 

-0.014 

0.050 

0.010 

0.010 

0.075 

0.084 

0.084 

0.100 

0.195 

0.195 

0.125 

0.325 

0.325 

0.150 

0.457 

0.457 

0.175 

0.574 

0.574 

0.200 

0.666 

0.666 

0.225 

0.725 

0.725 

0.250 

0.746 

0.746 

0.275 

0.725 

0.725 

0.300 

0.666 

0.666 

0.325 

0.574 

0.574 

0.350 

0.457 

0.457 

0.375 

0.325 

0.325 

0.400 

0.195 

0.195 

0.425 

0.084 

0.084 

0.450 

0.010 

0.010 

0.475 

-0.014 

-0.014 

0.500 

0. 

0. 
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The  base  problem  dependence  manifests  itself  in  the  w and  v matrices.  The 
inner  summation  of  these  matrices  can  go  no  higher  than  the  number  of  base 
problem  constants.  Convergence  is  determined  as  before.  Additionally,  and  i/J 
must  equal  w*  and  v\  respectively.  It  was  found  that  75  base  constants  were  required 
to  converge  even  a 10x10  v\,  at  which  point  the  10,10  element  still  differed  from  the 
melt  by  4.5%.  This  poor  convergence  is  directly  attributable  to  the  terms. 

To  improve  this  convergence,  a Shanks’  approximation  was  used14.  It  was  found  that 
15  terms  led  to  a convergence  better  than  that  of  90  terms  without  the 
approximation.  However,  the  exact  improvement  available  in  using  this 
approximation  is  indeterminable.  This  approximation  enabled  a matrix  size  of  26x10 
to  be  converged.  These  converged  elements  of  the  w and  v matrices  are  themselves 
summed  in  the  e matrix.  In  fact,  the  e matrix  must  be  doubly  summed.  The  inner 
sum  of  e was  never  satisfactorily  converged.  This  was  due  to  the  nature  of  the 
matrices.  The  w1  were  strictly  diagonally  dominate  and  the  diagonal  elements  were 
negative  and  the  off  diagonal  elements  were  positive.  Therefore,  the  inner  sum  of 
e is  positive  until  the  diagonal  of  is  encountered,  at  which  point  the  sum  becomes 
negative  and  then  once  again  becomes  more  positive  as  additional  terms  are  added. 
This  behavior  leads  to  a very  slowly  convergent  series  and  is  not  compatible  with  the 
Shanks’  approximation.  The  use  of  90  base  constants  with  the  Shanks’  approximation 
in  determining  matrices  still  failed  to  converge  this  inner  summation  of  e.  It 
seems,  therefore,  that  only  the  brute  force  method,  with  its  resulting  high  cost  of 
computer  time,  of  ever  increasing  the  number  of  base  constants  must  be  used.  In  an 
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effort  to  avoid  this,  we  return  to  the  defining  equations  to  investigate  the  source  of 
the  poor  convergence. 

Doing  this  one  finds  that  the  lack  of  convergence  is  attributable  to  terms 
which  are  associated  with  the  z-functions  of  Y.  These  z-functions  directly  arise  from 
the  no-slip  boundary  conditions  at  z-  0,0.  No-slip  is,  however,  only  one  of  the 
limiting  cases  of  closed  flow  problems.  The  no-stress  boundary  condition  at  z = 0,0 
also  results  in  closed  flow.  Thus,  the  no-stress  case  should  still  provide  the  desired 
insight  into  the  problem.  It  has  the  distinct  advantage  of  not  giving  rise  to  the  z- 
functions  found  to  be  so  troublesome  in  the  no-slip  case  and  should  not  arise  in  the 
convergence  problems  associated  with  no-slip.  It  is  disadvantageous  in  that  it  is  not 
as  physically  realistic  and  may  result  in  predicting  oscillations  where  none  are  found 
in  an  experimental  system.  Keeping  this  in  mind,  we  now  examine  the  no-stress 
system. 

The  equations  given  for  the  no-slip  system  are  directly  convertible  into 
equations  applicable  in  the  no-stress  system  by  allowing  the  F 3 and  F4  terms,  and  any 
further  terms  arising  from  them,  to  equal  zero.  This  results  in  a vastly  simplified  set 
of  equations  in  which  there  is  no  matrix  inversion  and  only  limited  summation. 
Specifically,  the  base  problem  reduces  to 

OO 

X nm  _ r>  nmT7mZm 

ts  1 n 
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where  the  constants  are  now  exactly  determined. 
The  perturbed  problem  reduces  to 
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Equations  (4.6.1),  (4.6.4),  and  (4.6.5)  are  still  applicable.  The  alphas  are  now  given 
by 
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a2 kl 
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and  the  equation  for  a™kl  does  not  change.  The  w’ s,  v’s,  and  e are  given  by 
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Equations  (4.7)  and  (4.8)  still  define  the  elements  of  A and  B.  The  a’s  change  as 
follows 
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where  y is  unchanged  but 
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and  B remains  unchanged.  Equations  (5.4)-(5.29)  define  the  no-stress  model  and  all 
previously  discussed  points  are  applicable  to  it. 

As  expected,  convergence  did  not  pose  a limitation  in  the  no-stress  case. 
Summations  only  occur  in  calculating  the  w’s,  v’s,  y’s,  and  in  the  elements  of  the 
eigenmatrix.  The  w,  v,  and  y summations  converged  by  the  25th  term  for  all 
parameter  sets.  The  eigenmatrix  proved  more  difficult  to  converge,  requiring  a 
summation  of  200  terms.  Tables  5.6-5.12  show  the  convergence  of  the  eigenvalue  as 
the  order  of  the  eigenmatrix  was  increased.  For  each  parameter  set  studied,  this 
convergence  was  rapid  and  an  order  of  20x20  was  deemed  sufficient  for  an  accurate 
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determination  of  the  eigenvalue.  The  boundary  conditions,  for  a typical  set  of 
parameters,  are  given  in  Tables  5.13-5.16.  The  boundary  conditions  were  well 
satisfied  except  for  the  normal  component  of  momentum,  which,  as  expected,  was 
dismally  poor. 

Tables  5.17-5.19  indicate  how  Ca”1  changes  while  varying  s,  0,  and  r j 
respectively.  It  is  seen  that  Ca”1  decreases  as  s decreases  for  an  aspect  ratio  which 
is  below  the  Rayleigh  limit.  Recalling  the  previous  reasoning  concerning  the  region 
of  stability  indicates  that  stability  is  increased  by  decreasing  s.  The  indicated  values 
of  the  critical  surface  tension,  ctc,  were  determined  from 


ac  = Ca”1  iim  vwall 


= Ca"1  (2.79mPa-s) 


3.58 


|im 


= CaT'-lO-5 

cm 


where  the  indicated  value  of  nm  is  for  GaAs  at  its  melting  point15.  The  maximum 
value  of  5 = 1500  was  chosen  because  this  is  approximately  the  ratio  of  viscosities  of 
B203  to  GaAs16.  It  is  also  determined  that  increasing  the  melt  height  for  a fixed 
base  problem  interface  location  and  viscosity  ratio  results  in  an  increase  in  stability 
whereas  increasing  the  base  interface  location  for  a fixed  height  and  viscosity  ratio 
results  in  a decrease  in  stability.  The  imaginary  result  for  0 = 1.  indicates  that  the 
system  is  oscillatory  for  all  surface  tensions. 
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TABLE  5.6  Eigenvalue  convergence  for  rx  = .5,  0=3.,  and  5 = 1500 


Eigenmatrix  Order 

Ca1 

2x2 

.9018 -105 

4x4 

.6603 -105 

6x6 

.6495 -105 

8x8 

.6466 -105 

10  x 10 

.6459 -105 

12  x 12 

.6453 -105 

14  x 14 

.645 1-105 

16  x 16 

.6450-105 

18  x 18 

.6450-105 

20x20 

.6449 -105 

TABLE  5.7  Eigenvalue  convergence  for  ^ = .5,  0=3.,  and  5=150 


Eigenmatrix  Order 

Ca1 

2x2 

.9104-104 

4x4 

.6675 -104 

6x6 

.6567 -104 

8x8 

.6538 -104 

10  x 10 

.6530 -104 

12  x 12 

.6524 -104 

14  x 14 

.6523 -104 

16  x 16 

.6522 -104 

18  x 18 

.6521 -104 

20x20 

.6520 -104 
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TABLE  5.8  Eigenvalue  convergence  for  rx  = .5,  (3=3.,  and  5 = 15. 


Eigenmatrix  Order 

Ca1 

2x2 

.1008  TO4 

4x4 

.7508 -103 

6x6 

.7388 -103 

8x8 

.7359 -103 

10  x 10 

.7350-103 

12  x 12 

.7345 -103 

14  x 14 

.7343 -103 

16  x 16 

.7342 -103 

18  x 18 

.7341 -103 

20  x 20 

.7341 -103 

TABLE  5.9  Eigenvalue  convergence  for  ^ = .5,  /3=3.,  and  5 = 1.5 


Eigenmatrix  Order 

Ca1 

2x2 

.3036 -103 

4x4 

.2546 -103 

6x6 

.2530T03 

8x8 

.2527 -103 

10  x 10 

.2526 -103 

12  x 12 

.2525 -103 

14  x 14 

.2525 -103 

16  x 16 

.2525 -103 

18  x 18 

.2525 -103 

20  x 20 

.2525  TO3 
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TABLE  5.10  Eigenvalue  convergence  for  ^ = .5,  13=2.,  and  5=1500. 


Eigenmatrix  Order 

Ca1 

2x2 

.1935  TO5 

4x4 

.1806-105 

6x6 

.1800  TO5 

8x8 

.1800-105 

10  x 10 

.1799-105 

12  x 12 

.1799T05 

14  x 14 

.1799-105 

16  x 16 

.1799-105 

18  x 18 

.1799T05 

20  x 20 

.1799-105 

TABLE  5.11  Eigenvalue  convergence  for  ^ = .6,  f3=  3.,  and  5 = 1500. 


Eigenmatrix  Order 

Ca1 

2x2 

.1760-107 

4x4 

.1632-107 

6x6 

.1628-107 

8x8 

.1627-107 

10  x 10 

.1626  TO7 

12  x 12 

.1626  TO7 

14  x 14 

.1626  TO7 

16  x 16 

.1626  TO7 

18  x 18 

.1626  TO7 

20  x 20 

.1626  TO7 
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TABLE  5.12  Eigenvalue  convergence  for  ^ = .55,  (3=3 .,  and  s=  1500. 


Eigenmatrix  Order 

Ca1 

2x2 

.1743  TO6 

4x4 

.1517-106 

6x6 

.1509  TO6 

8x8 

.1507  TO6 

10  x 10 

.1505  TO6 

12  x 12 

.1505 -106 

14  x 14 

.1505 -106 

16  x 16 

.1505  TO6 

18  x 18 

.1505  TO6 

20x20 

.1504  TO6 
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TABLE  5.13  Moving  wall  no-slip  boundary  condition  for  ^ = .5,  /?  = 3.,  5=  1500.  . 


z 

& 

Ti 

1 

r=\ 

dr 

0.00 

0. 

0.15 

0.28  TO'13 

0.30 

0.56 -10'13 

0.45 

0.81-10"13 

0.60 

1.04 -10' 13 

0.75 

1.23-10”13 

0.90 

1.39-10'13 

1.05 

1.5  MO'13 

1.20 

1.59  TO' 13 

2.25 

1.65-10"13 

1.35 

rH 

1 

o 

r— < 

r- 

\o 

1.5 

1.66  TO'13 

1.65 

1.61  TO'13 

1.80 

1.52-10"13 

1.95 

1.39T0'13 

2.10 

1.23 -10’ 13 

2.25 

1.03 -10'13 

2.40 

0.81  TO'13 

2.55 

0.55 -10- 13 

2.70 

0.28 -10' 13 

2.85 

0.28 -10' 13 

3.00 

0. 
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TABLE  5.14  Interfacial  no-slip  boundary  condition  for  ^ = .5,  P=3.,  5 = 1500.  . 


z 

m 

*(0)  + 

dT(0) 

r=rl 

^(0) 

8(0) 

r=ri 

dr 

dr 

0.00 

0. 

0. 

0.15 

0.253 

0.252 

0.30 

-0.562 

-0.563 

0.45 

-1.135 

-1.135 

0.60 

-1.483 

-1.483 

0.75 

-1.707 

-1.707 

0.90 

-1.835 

-1.835 

1.05 

-1.872 

-1.872 

1.20 

-1.826 

-1.825 

1.35 

-1.706 

-1.706 

1.50 

-1.532 

-1.532 

1.65 

-1.326 

-1.326 

1.80 

-1.107 

-1.107 

1.95 

-0.889 

-0.888 

2.10 

-0.688 

-0.688 

2.25 

-0.517 

-0.516 

2.40 

-0.378 

-0.378 

2.55 

-0.266 

-0.266 

2.70 

-0.180 

-0.181 

2.85 

-0.080 

-0.081 

3.00 

0. 

0. 
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TABLE  5.15  Interfacial  tangential  stress  boundary  condition  for  ^ = .5,  /J  = 3., 
s=  1500.. 


z 

. m 1 d^m 

h a + + v ' 

ri 

s 

S2T'i13T'/J([)  < 1 

\ (0)  dr2  r a- 

Jrl 

0.0 

0. 

0. 

0.15 

-9.54 

-9.49 

0.30 

11.70 

11.73 

0.45 

20.69 

20.71 

0.60 

25.12 

25.13 

0.75 

27.81 

27.81 

0.90 

29.26 

29.25 

1.05 

29.15 

29.14 

1.20 

27.98 

27.97 

1.35 

25.37 

25.35 

1.50 

21.91 

21.89 

1.65 

18.08 

18.07 

1.80 

14.42 

14.40 

1.95 

10.63 

10.61 

2.10 

7.48 

7.47 

2.25 

5.04 

5.04 

2.40 

3.54 

3.54 

2.55 

2.32 

2.33 

2.70 

1.81 

1.83 

2.85 

0.39 

0.42 

3.00 

0. 

0. 
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TABLE  5.16  Interfacial  normal  stress  boundary  condition  for  ^ = .5,  p =3.,  s=  1500. 


z 

Ihs  of  equation  (3.7.2) 

rhs  of  equation  (3.7.2) 

0.00 

-3.66T0+4 

-0.0001  T0+5 

0.15 

0.42 -10 + 4 

-0.40  TO + 5 

0.30 

0.13-10+4 

-0.14T0+5 

0.45 

-0.20-10+4 

-4.24  TO + 5 

0.60 

-0.37  TO + 4 

0.14-10+5 

0.75 

-0.48  TO + 4 

0.19T0+5 

0.90 

-0.57 -10 + 4 

0.30-10+5 

1.05 

-0.63 -10 + 4 

0.37-10+5 

1.20 

-0.67T0+4 

0.56  TO + 5 

1.35 

-0.68T0+4 

0.62T0+5 

1.50 

-0.66T0+4 

0.65-10+5 

1.65 

-0.61  TO + 4 

0.62-10+5 

1.80 

-0.55  TO + 4 

0.63  TO + 5 

1.95 

-0.48  TO + 4 

0.46T0+5 

2.10 

-0.39  TO + 4 

0.35 -10+5 

2.25 

-0.31T0+4 

0.23  TO + 5 

2.40 

-0.24  TO + 4 

0.24  TO + 5 

2.55 

-0.19T0+4 

0.13T0+5 

2.70 

-0.14-10+4 

5.33*10+5 

2.85 

-0.14T0+4 

-0.41  T0+5 

3.00 

0.13T0+4 

0.0002  TO + 5 
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TABLE  5.17  The  affect  of  varying  s on  ac 


rl 

0 

s 

Ca^1 

CTc 

(dyne/cm) 

.5 

3. 

1500. 

64490. 

.64 

.5 

3. 

150. 

6520. 

.065 

.5 

3. 

15. 

734.1 

.0073 

.5 

3. 

1.5 

252.5 

.0025 

TABLE  5.18  The  affect  of  varying  (3  on  ac 


rl 

0 

s 

Ca;1 

ac 

(dyne/cm) 

.5 

3. 

1500. 

64490. 

.64 

.5 

2. 

1500. 

17990. 

.18 

.5 

1. 

1500. 

Im 

Im 

TABLE  5.19  The  affect  of  varying  r1  on  ac 


ri 

0 

s 

Ca-1 

ac 

(dyne/cm) 

.5 

3. 

1500. 

64490. 

.64 

.55 

3. 

1500. 

150400. 

1.5 

.6 

3. 

1500. 

1626000. 

16.3 
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Tables  5.20-5.22  indicate  the  levels  of  base  problem  normal  stress  at  the 
interface  for  the  comparisons  made  in  Tables  5.17-5.19.  From  Table  5.20  it  is  seen 
that  the  level  of  stress  increases  for  a decrease  in  5.  The  fact  that  this  also  results 
in  an  increase  in  stability,  as  indicated  by  the  decreasing  Ca~  \ is  not  fully  understood 
but  seems  to  indicate  that  both  the  interface  deflection  and  the  perturbed  problem 
normal  stress  play  important  roles  in  determining  stability.  Table  5.21  again  indicates 
that  an  increase  in  base  normal  stress  results  in  stability.  Table  5.22,  however,  leads 
one  to  the  opposite  conclusion.  In  this  case  it  is  felt  that  the  increasing  aspect  ratio 
of  the  encapsulant  is  having  an  effect  on  the  stability  of  the  interface.  The  aspect 
ratio  of  the  encapsulant,  heretofore  ignored,  does  seem  to  play  a critical  role.  Table 
5.23  shows  the  base  problem  interfacial  no-slip  boundary  condition  for  ^ = .7,  /?=  3., 
and  s = 1500.  which  is  an  encapsulant  aspect  ratio  of  10.  Physical  insight  suggests,  and 
the  values  given  in  Table  5.4  concur,  that  the  value  of  v,  should  be  negative  for  wall 
movement  in  the  +z  direction.  However,  Table  5.23  indicates  positive  and  negative 
vz  f°r  positive  wall  movement.  This  may  only  occur  if  more  than  one  vortex  cell, 
such  as  depicted  in  Figure  5.2,  arises  in  the  encapsulant.  For  such  a cell  pattern  to 
occur,  no-slip  must  be  violated  somewhere.  It  appears  therefore,  that  the  equations 
can  no  longer  depict  an  accurate  model  of  the  system  for  these  large  aspect  ratios. 
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TABLE  5.20  Base  Problem  Normal  Stress  for  varying  viscosity  ratios 


z 

/ 

d1wm  d2?*  , . , 

value  of  - s for  s equal  to 

drdz  drdz  . 

T\ 

1500. 

150. 

15. 

1.5 

.3 

0.688 

0.694 

0.698 

0.728 

.6 

0.198 

0.206 

0.225 

0.358 

.9 

0.089 

0.096 

0.109 

0.208 

1.2 

0.035 

-0.042 

0.049 

0.099 

1.5 

-0.009 

-0.001 

0.000 

0.000 

1.8 

-0.053 

-0.044 

-0.049 

-0.099 

2.1 

-0.110 

-0.098 

-0.110 

-0.208 

2.4 

-0.228 

-0.209 

-0.225 

-0.358 

2.7 

-0.745 

-0.699 

-0.698 

-0.728 
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TABLE  5.21  Base  Problem  Normal  Stress  for  varying  zone  heights 


z 

\ 

value  of  - s for  0 equal  to 

drdz  drdz  _ M 

rl 

3. 

2. 

1. 

.1 

- 

- 

0.120 

.2 

- 

0.991 

0.856 

.3 

0.688 

- 

0.494 

.4 

- 

0.396 

0.205 

.5 

- 

- 

-0.026 

.6 

0.198 

0.147 

-0.260 

.7 

- 

- 

-0.559 

.8 

- 

0.050 

-0.945 

.9 

0.089 

- 

-0.137 

1. 

- 

-0.013 

- 

1.2 

0.035 

-0.078 

- 

1.4 

- 

-0.180 

- 

1.5 

-0.009 

- 

- 

1.6 

- 

-0.441 

- 

1.8 

-0.053 

-0.108 

- 

2.1 

-0.110 

- 

- 

2.4 

-0.228 

- 

- 

2.7 

-0.745 

- 

- 
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TABLE  5.22  Base  Problem  Normal  Stress  for  varying  interface  locations 


z 

' _ ' 

a2?"1 

value  of  - s for  r,  equal  to 

drdz  drdz  J 14 

r\ 

0.50 

0.55 

0.60 

.3 

0.688 

0.994 

3.312 

.6 

0.198 

0.422 

2.335 

.9 

0.089 

0.241 

1.592 

1.2 

0.035 

-0.109 

0.806 

1.5 

-0.009 

-0.009 

-0.009 

1.8 

-0.053 

-0.128 

-0.824 

2.1 

-0.110 

-0.262 

-1.614 

2.4 

-0.228 

-0.451 

-2.364 

2.7 

-0.745 

-0.105 

-3.369 
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TABLE  5.23  Interfacial  no-slip  boundary  condition  for  ^ = .7,  (3  =3.,  and  5 = 1500.  . 


z 

dWm 

* Jrx 

dWe 

<>  Jr, 

0.0 

0. 

0. 

0.3 

-0.7290 

-0.7290 

0.6 

-0.5369 

-0.5369 

0.9 

0.1511 

0.1511 

1.2 

0.8833 

0.8833 

1.5 

1.193 

1.193 

1.8 

0.8833 

0.8833 

2.1 

0.1511 

0.1511 

2.4 

-0.5369 

-0.5369 

2.7 

-0.7290 

-0.7290 

3.0 

0. 

0. 
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FIGURE  5.2  Multiple  cell  pattern  in  which  no-slip  is  violated  for  the  lower  cell 


CHAPTER  6 

CONCLUSIONS  AND  SCOPE 


The  stability  of  liquid  encapsulated  melt  zones  (LEMZ)  is  of  interest  to  a 
wide  range  of  engineers,  but  particularly  those  concerned  with  crystal  growth.  In  this 
study  analytical  solutions  were  developed  for  a model  problem  LEMZ  under  both 
no-slip  and  no-stress  boundary  conditions.  The  no-stress  solutions  were  evaluated 
and  critical  inverse  capillary  numbers  were  determined.  From  this  the  value  of  the 
critical  surface  tension  can  be  determined  for  a given  encapsulant  and  melt  fluid.  It 
is  dependent  upon  the  three  parameters  of  the  model  rv  the  location  of  the  base 
problem  vertical  interface,  /3,  the  height  of  the  zone,  and  s,  the  ratio  of  encapsulant 
to  melt  viscosities. 

It  was  determined  that  for  aspect  ratios  below  the  Rayleigh  limit  that  the 
stability  of  the  zone  decreases  as  5 is  increased.  This  was  the  major  objective  of  the 
study.  At  the  same  time,  it  was  found  that  the  level  of  normal  stress  due  to  the  base 
problem  flow  also  was  decreasing.  This  leads  one  to  the  surprising  conclusion  that 
an  increase  in  stress  increases  the  stability.  Further  study  of  this  phenomena,  to 
determine  the  affect  of  the  perturbed  normal  stress  and  the  deflection  of  the 
interface,  are  called  for  before  it  can  be  explained.  It  should  be  stressed  that  this 
would  only  help  in  predicting  the  trend  of  the  results.  A calculation  would  still  be 
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required  to  determine  oc.  Extending  the  analysis  to  aspect  ratios  above  the  Rayleigh 
limit  is  also  called  for. 

Secondary  points  of  interest  were  how  changes  in  r1  and  /9  affect  stability.  It 
was  found  that  decreasing  /3  for  fixed  and  s increased  stability  whereas  increasing 
ri  for  fixed  /3  and  s decreased  stability.  This  is  interpreted  as  meaning  lower  melt 
aspect  ratios  result  in  greater  stability.  This  also  requires  further  study,  as  it  seems 
that  how  the  aspect  ratio  is  obtained,  either  by  changing  r ^ or  /3,  can  also  affect  the 
results.  This  led  to  one  of  the  most  surprising  phenomena  in  this  study,  the  inability 
of  the  model  to  accurately  depict  the  system  for  high  encapsulant  aspect  ratios. 
Exactly  what  is  occurring  when  this  happens  needs  to  be  understood,  although  it 
appears  that  the  flow  is  unable  to  turn  for  these  high  ratios. 

It  is  important  to  remember  that  these  results  are  only  applicable  in  the 
region  of  steady  solutions  generated  in  this  study.  No  comparison  can  be  made 
between  the  ac  determined  here  and  the  oc  which  would  result  if  time  dependency 
is  allowed  ( in  the  kinematic  condition).  This  steady  state  assumption  can  be  easily 
relaxed  in  further  work.  The  other  assumptions  of  the  model  were:  zero  gravity, 
isothermal,  creeping  flow,  and  axi-symmetric  disturbances.  Each  of  these  can  also 
be  relaxed.  The  usefulness  of  relaxing  creeping  flow  is  not  apparent  as  this 
accurately  depicts  crystal  growth.  Relaxing  axi-symmetric  disturbances  forces  one  to 
a purely  numerical  solution  to  solve  the  resulting  three  dimensional  flow  and  this  is 
undesirable.  Relaxing  the  zero  gravity  assumption,  either  by  allowing  a non-zero 
value  of  gravity  or  by  allowing  a density  difference  between  the  fluids  in  an 
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experiment,  is  useful  if  one  is  interested  in  earth  processing  but  is  not  as  necessary 
if  the  primary  interest  is  space  processing.  The  isothermal  assumption  does  deserve 
to  be  relaxed  in  future  work.  Doing  so  will  introduce  another  convective  flow,  the 
Marangoni  effect.  It  is  likely  that  this  will  affect  the  LEMZ  stability  but  it  is  unclear 
how  it  will  affect  the  conclusions  with  regards  to  s. 

Theoreticians  must  constantly  be  reminded  that  theory  demands  experimental 
verification.  It  is  a void  in  this  study  which  needs  to  be  filled.  One  can  envision  a 
possible  experimental  apparatus  along  the  lines  of  a Plateau  chamber.  It  would 
require  modification  to  allow  for  wall  movement  and  equal  encapsulant  and  melt 
heights.  Water  and  silicone  oil  may  prove  to  be  suitable  fluids  for  such  an 
experiment.  They  are  immiscible,  the  density  of  the  silicone  oil  can  be  adjusted  to 
achieve  a neutral  buoyancy  affect  to  simulate  zero  gravity,  and  silicone  oil  is  available 
in  different  viscosities,  one  of  which  results  in  a silicone  oil  to  water  ratio  of 
approximately  the  same  as  the  GaAs  to  B203  ratio  of  1500.  . 

Finally,  the  no-slip  case  deserves  to  be  revisited.  Although  it  is  not 
anticipated  to  provide  much  more  insight  as  to  the  role  of  s,  it  is  the  more  physically 
realistic  case.  From  a theoretical  point  of  view,  a deeper  understanding  of  the 
convergence  problems  may  help  one  to  avoid  them  in  the  future. 


APPENDIX  A 

DERIVATION  OF  DIMENSIONLESS  EQUATIONS 


Domain  Equations 

The  equation  of  motion  for  an  incompressible  fluid  in  the  absence  of  gravity 


is 


— = -Iw  - !v-r  - v-w  (a.i) 

ar  p p 


Applying  Newton’s  law  of  viscosity  results  in 


dv 

dt 


-lw  + vV 
P 


Vv  + (W^  - v • W 


(A.2) 


Introducing  the  dimensionless  variables 


yields 


vc  dv  * 
dt* 


(A.3) 


or 
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where 


Strouhal  number  = Sr  = — L_ 

Vc 


Reynolds  number  = Re  = 


rcvc 


Choosing  P c = 


li  rv 


C'  C 


yields 


Re  Sr 


dv* 


dt* 


Rev*  • V*v* 


(A.5) 


and  assuming  steady  state  and  creeping  flow  results  in 


V*P* 


V*v 


M*V* 


(A.6) 


Equation  (A.6)  could  also  be  obtained  for  the  unsteady  problem  in  the  creeping  flow 
regime  by  choosing  tc  so  that  Sr  = 1 and  assuming  the  left-hand-side  of  equation  (A.5) 
to  equal  zero  due  to  the  small  Reynolds  number.  Since  all  quantities  are  now 
dimensionless  the  will  be  dropped,  but  understood,  from  this  point  onward. 
Equation  (A.6)  is  the  familiar  Navier-Stokes  equation  for  an  incompressible  fluid 
with  constant  properties  in  dimensionless  form. 
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Due  to  the  absence  of  gravitational  effects  and  the  cylindrical  geometry,  it  is 
possible  to  assume  that  v&  and  partial  derivatives  with  respect  to  0 equal  zero. 
Therefore,  only  the  r and  z components  of  equation  (A.6)  remain.  These  are 


dP 


/ 


1 

r dr 


(A.7.1) 


dP 

Iz 


/ 


\ 


1 & Z x &vz 

r to  ~a? 


Since 


d1P  = &P 
dzdr  drdz 


(A.7.2) 


(A.8) 


differentiating  equation  (A.7.1)  with  respect  to  z and  equation  (A.7.2)  with  respect 
to  r and  equating  the  resultant  equations  yields 


d 

/ \ 
d*vr  1 dvr  vr  d2vr 

d 

' \ 

d^z  1 toz  02vz 

& 

dr2  r dr  f2  &2 

dr 

+ + 

,dz2  r to  0,2  ^ 

(A.9) 


or 


dr* 


dz  dr 


+ \d_ 
r dr 


to r _ 
dz  dr 


dvr  dvz 
dz  dr 


dz2 


dvr  dvz 
dz  dr 


= o (A.10) 


Defining 
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results  in 


0)  = 


for  _ & z 
dz  dr 


^co  i 3o)  _ <i>  a2* 

dr2  r fo  ~r2+~S2 


= 0 


or 


where 


L(rw)  = 0 


L 3 _ J.  A.  + ^L 


,2  r 6lr 


Introducing  the  stream  function,  Y,  such  that 


1 dW  , 1 dl? 

v,  = and  v.  = 


r dr 


r dz 


(thereby  satisfying  continuity)  and  rewriting  <0  in  terms  of  Y yields 


0)  = 


l d2?  _ 1 8Y  Id2!! 
r dr2  r2  fo  + r dz2 


(A.  11) 


(A.  12) 


(A.  13) 


Substituting  this  into  equation  (A.  12)  results  in  L2Y  = 0 or 
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+ 2_^Z_  + 2 2 ^ + 3 ^ 3 dW 

dr 4 dr^dz^  dz ^ r dr ^ r drdz^  r ^ ^r 


(A.  14) 


and  it  follows  that  LY  = rto.  This  describes  the  variation  of  Y with  respect  to  r and 
z.  It  is  a fourth  order  partial  differential  equation  in  terms  of  r and  z.  Equation 
(A.  14)  is  equally  applicable  in  the  both  the  melt  and  the  encapsulant  phases,  with  the 
only  difference  being  the  boundary  conditions  which  are  applied.  In  each  liquid  eight 
boundary  conditions  are  required. 

Interfacial  Equations 

The  equations  of  interest  are  the  momentum  balance  and  the  kinematic 
condition. 

Momentum  Balance. 

The  interfacial  momentum  balance  for  incompressible  fluids  is 


i m 


n 


- nmn 


W 


m 


= Pen 


lien 


- 2Hon 


(A.  15) 


where  n is  the  normal  interfacial  surface  vector  and  H is  the  mean  curvature  of  the 
interface  (See  Appendix  B).  Using  dimensionless  variables  as  before  results  in 


Pmn 


n 


W 


m 


sn  • 


2//Ca-1n 


(A.  16) 


where  the  * has  been  dropped  and  Capillary  number  = Ca  = 


Hmv 


c and  s=- — 


,m 
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Kinematic  Condition. 

The  kinematic  condition  describes  the  velocity  of  the  interface.  It  is  by 
definition 


To  eliminate  the  surface  velocity,  i/,  from  equation  (A.  17),  the  interfacial  mass 
balance  is  used,  with  the  restriction  of  no  mass  transfer  between  the  two  liquids. 
This  results  in 


Therefore,  v/n=us=tf.  Using  dimensionless  variables  as  before  results  in 


(A.  18) 


The  no  slip  conditions  dictate  that 


r 


v 


r 


m 


(A.  19.1) 


Dr? 


v 


e 


(A19.2) 


Dr 


r 


r 


where  once  again  the  * has  been  dropped. 


APPENDIX  B 

DETERMINATION  OF  SURFACE  VECTORS  AND  MEAN  CURVATURE 


Surface  Vectors 

Normal  Vector 

In  general,  a surface  in  cylindrical  coordinates,  is  described  by  r?=f(0^)  or 
g(rj,Q1z)  = ri-i.  The  normal  to  this  surface  is  given  by 

n - (B.1) 

Pg  ■ Vg 


which  evaluates  to 


n = 

i ♦ fiiif  ♦ 

'<V 

2 

_ 1 
2 

1 dr]  dr] 

e r -ea  - — ~e. 

U 3e  J 

/ T]  dd  e dz  z 

(B.2) 


If  axial  symmetry  is  assumed  equation  (B.2)  reduces  to 


n = 

1 + 

fan] 

2 

_ 1 
2 

dr] 

er  - — ~ey 

l&J 

[ & z\ 

(B.3) 


It  should  be  noted  that  equation  (B.3)  is  valid  in  dimensionless  form  as  well. 
Tangent  Vectors 

The  two  tangent  vectors  to  the  surface  must  be  determined  so  that 
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(B.4.1) 

(B.4.2) 

(B.4.3) 


n-tl  = 0 
n -t2  = 0 

h 'h  = 0 


are  satisfied.  Since  t1  and  t2  each  have  3 components,  there  are  6 unknowns  in 
equations  (B.4).  Therefore,  3 of  the  unknowns  are  free  to  be  specified.  The 
resultant  tangent  vectors  depend  upon  the  choice  of  these  unknowns.  Writing  tl  and 

h as 

h = her  + fad  + l\ez 

h = *2er  + llee  + (2ez 


and  substituting  into  equations  (B.4)  results  in 


l 

l 


2 1 dp 


(B.5.1) 


t 


1 

2 


2 1 dp 

2 p de 


,3  dp 
2^ 


= 0 


(B.5.2) 


t2tl 

hl2 


.2.2 
+ *1*2 


3 3 

hh 


= o 


(B.5.3) 


Choosing 


1 drl  ,2  _ _dp  f3  _ drl 
1 = ■¥  ’ 1 = "■&  ’ '2  = &■ 


substituting  into  equations  (B.5)  and  solving  for  the  unknowns  yields 
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'l  - 


t'l  ~ 


l ♦ ±^L 

rj  50 

1 - 1 dr} 


1 - 


rj  50 

1 5rj 

~nle 


-i 


i 1 3rj 
1 + L + 

r)  50 


(5r^2 

dz 


1 dr]  ( 1 dr]  Y2 


r?  50 


r?  50 


ranf 
* , 


Therefore 


l 


J_5rf 

n ae. 


(B.6) 


and 


h ~ 


1 - 


1 dr] 

-1 

1 dr] 

1 + 

' 1 dr]' 

r]  50 

r]  50 

30, 

5r; 


1_  1 

-1 

i 1 3r? 

1 + L + 

'drT 

2 

n 30 

rj  50 

c 0 

dr] 

+ — e7 
dz  z 


(B.7) 


It  is  easily  verified  that  equations  (B.4)  are  satisfied.  The  choice  of  t\,  t\,  and  t\  is 
made  so  that  if  r]=rl  then  n=er  ti=ez,  and  t2=e 6,  which  corresponds  to  the 
interfacial  surface  vectors  of  the  base  problem  interface.  If  axial  symmetry  is 
assumed  then  equations  (B.6)  and  (B.7)  reduce  to 


+ e. 


(B.8) 
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+ 


(B.9) 


It  is  noted  that  equations  (B.8)  and  (B.9)  are  valid  in  dimensionless  form. 

Mean  Curvature 

The  mean  curvature  of  an  interface  is  given  by 

-2 H = Vs-n  (B.10) 


where 


i a a 

—ea — + e7  — 

t)  e ae  z dz 


Using  equation  (B.2)  in  equation  (B.10)  results  in 


-2 H = 


7 

1 + 

'dn  V 

„2  * f*Lf 

_ 3 
2 ■ 

1 + 

'dq' 

2 

J l aej  J 

+ 


2 drj  dy  a2^ 

“ dz  ae  n dGdz 


dr!  dr! 

aF  ae 


(B.ll) 


+ 


f*j)2 . 

,<»JJ 


If  axial  symmetry  is  assumed  then  equation  (B.ll)  reduces  to 
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-2 H = 


1 + 


fan 


\2 


1 cPr) 


1 dz ' 


1 + 


dr? 

at 


\2 


dr?  )2  cPr, ? 


at 


ct' 


It  is  noted  that  equation  (B.12)  is  valid  in  dimensionless  form. 


(B.12) 


APPENDIX  C 

BOUNDARY  CONDITIONS 


Equation  (A.  14)  describes  the  flow  for  both  the  melt  and  the  encapsulant.  It 
is  a fourth  order  partial  differential  equation  in  two  variables.  Therefore,  8 boundary 
conditions,  4 for  each  variable,  are  required  in  both  the  melt  and  the  encapsulant. 
Additionally,  the  location  of  the  liquid-liquid  interface  is  unknown  so  another 
boundary  condition  is  required,  for  a total  of  17. 

Of  these  17,  10  are  the  usual  conditions  of  no  slip  and  no  flow  at  the  solid 
surfaces.  At  the  centerline,  the  conditions  of  no  flow  and  bounded  velocity  exist. 
These  12  boundary  conditions  are  summarized  in  Table  C.l.  The  5 remaining 
conditions  are  specified  at  the  liquid-liquid  interface.  These  are  the  kinematic 
condition  of  the  melt,  the  kinematic  condition  of  the  encapsulant,  no  slip,  the 
continuity  of  tangential  stress  and  the  normal  component  of  the  interfacial 
momentum  balance. 
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TABLE  C.l  Non-Interfacial  Boundary  Conditions 


MELT 

ENCAPSULANT 

vf  (r,  z=0)  = 0 
\ (r,  z=0)  = 0 

v/  (r,  z=0)  = 0 

v/  (r,  z=0)  = 0 

v”  (r,  z=p)  = 0 
vrm  (r,  z=/9)  = 0 

vz  (r,  z=/3)  = 0 
v*  (r,  z=/3)  = 0 

vrm  (r=0,  z)  = 0 
\ 

dvf 

z = 0 

, & )r= 0 

v,*  (r=l,  z)  = 0 
v/  (r=l,  z)  = 1 

Interfacial  Boundary  Conditions 

Kinematic  Conditions 

Equation  (A.  19.1)  describes  the  kinematic  condition  for  the  melt.  Assuming 
steady  state  yields 


dr? 

dz 


Introducing  the  stream  function,  Y,  yields 


dr?  dWm  _ dWm 
dT  dr  dz 


(C.1) 


(C.2) 


The  kinematic  condition  for  the  encapsulant  is  analogous  to  this  and  is  given  by 

_ dp  dT  = dT^  (C.3) 

dz  dr  dz 


99 


No  Slip 


The  no  slip  condition  at  the  interface  dictates  that 


vm  = ve  -tx 


(C.3.1) 


vm  t2  = ve  ’t2 


(C.3.2) 


Using  equation  (B.8)  in  equation  (C.3.1)  and  equation  (B.9)  in  equation  (C.3.2) 
(thereby  assuming  axial  symmetry)  results  in 


dr]  m m dr)  e e 
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(C.4.2) 


It  is  apparent  that  equation  (C.4.2)  is  the  same  as  equation  (C.4.1),  so  there  is  only 
one  no  slip  condition.  Substituting  equation  (C.l)  and  its  equivalent  for  the 
encapsulant  into  equation  (C.4.1)  yields 


(C.5) 


Introducing  the  stream  function  results  in 

dVm  = (C.6) 

dr  dr 
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Tangential  Stress 

The  tangential  stress  is  obtained  by  taking  the  dot  product  of  the  interfacial 
momentum  balance  with  the  interfacial  surface  tangent  vectors.  Doing  this  for  the 
first  tangent  vector,  tv  yields 


n 


= sn  • 


(C.7) 


and  assuming  axial  symmetry  results  in 
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Using  t2  and  assuming  axial  symmetry  results  in 
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It  is  clear  that  when  axial  symmetry  is  assumed  equation  (C.9)  results  in  the  same 
condition  as  equation  (C.8).  Thus  there  is  only  one  tangential  stress  boundary 
condition,  as  was  the  case  for  no  slip. 

In  the  perturbed  problem,  the  assumption  of  axial  symmetry  is  fairly  severe. 
However,  as  it  is  valid  in  the  base  problem,  it  is  not  unreasonable  for  it  to  hold  in 
the  perturbed  problem.  In  order  to  relax  this  assumption,  the  full  3-dimensional 
problem  must  be  solved.  This  is  not  possible  analytically  and  would  require  a 
numerical  solution. 

Introducing  the  stream  function  into  equation  (C.8)  results  in 
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(C.10) 
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where 


^2 Tt^m  or  e 
frdz 


or  e 


has  been  used. 

Normal  Component 

The  normal  component  of  the  momentum  balance  is  obtained  by  taking  the 
dot  product  of  the  interfacial  momentum  balance  with  the  interfacial  surface  normal 


vector 
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Vvm  + (vvm)T\-n  = Pe  - sn 
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+ (we)T]  •/i-2//Ca“1  ^C11^ 
Assuming  axial  symmetry  and  using  equations  (B.3)  and  (B.12)  results  in 
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Rearranging  equation  (C.9)  to 
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and  substituting  this  into  equation  (C.12)  results  in 
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Introducing  the  stream  function  into  equation  (C.13)  results  in 
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has  again  been  used. 


APPENDIX  D 

ORTHOGONALITY  OF  P FUNCTIONS 


The  function  P^c^r),  defined  by 

PiM  - JiM  - (D-d 

Yi(Vi) 


can  be  shown  to  be  orthogonal  by  considering  the  characteristic  equation  (2.3).  One 
solution  to  this  equation  is  Pjfc,,/')  and  another  is  P^c^r).  From  equation  (2.3) 


Integrating  over  the  interval  to  1,  using  integration  by  parts  to  eliminate  the  second 
derivatives  and  simplifying  results  in 
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m 


~ cl)l1irPl[cmr)?l[cnr)dr  (D3) 


It  is  clear  that  Pi(cm/'1)  = P1(c/Jr1)  = 0 and  since  the  cm  and  cn  are  determined  such 
that  Pi(cm)  = P1(cn)  = 0 equation  (D.3)  reduces  to 


l\r?l[cmr)?lcnr)&  = 0 


(D.4) 


which  proves  the  orthogonality  of  the  functions  with  respect  to  the  weighting  function 


APPENDIX  E 

DETERMINATION  OF  THE  PERTURBED  PRESSURE  DIFFERENCE 

Equation  (A.6)  describes  the  velocity  in  both  fluid  phases.  Taking  the 
divergence  of  this  equation  results  in 

V-VP  = n V-V^  = n V^V-v) 
and  applying  continuity  yields 


Equation  (E.l)  is  valid  in  both  the  encapsulant  and  the  melt.  Perturbing  and 
equating  0(en)  terms  results  in 


The  perturbed  pressure  is  determined  by  solving  equation  (E.2)  subject  to  unknown 
Neumann  boundary  conditions.  Neumann  conditions  are  selected  so  that  equations 
(A.7)  may  be  used  to  introduce  the  stream  function. 


V2P  = 0 


(E.l) 


VZP  = 0 


(E.2) 
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How  these  boundary  conditions  are  formulated  is  critical  and  a detailed 
analysis  of  how  quantities  are  perturbed  and  where  they  are  evaluated  is  called  for. 
In  general,  a quantity  / (where  / may  be  a variable,  function,  or  functional)  is 
perturbed  in  the  following  manner 
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Specifically,  in  the  domain 
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where  the  interchangeability  of  the  order  of  differentiation  has  been  used.  Also 


V^rm  = V2!  — — 
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so 


(*H>] = ^ 


[0] 


(E.3) 


The  boundary  conditions  are  of  the  form 


dP[o\  , ap[0, 

= gl  and  -±1  = g2 


Introducing  the  stream  function  into  equations  (A.7)  results  in 
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r dr2  r2  * drdz: 


(E.4.2) 


Perturbing  the  left-hand-sides  of  equations  (E.4),  while  applying  the  same  argument 
as  was  used  for  V2P,  results  in 
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dP 

dz 
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+ € 

dz  [ & 

W + e_fP(o| 

dz  dz 


[0] 


Therefore,  the  O(e)  terms  of  the  perturbed  right-hand-sides  of  equations  (E.4)  define 
the  unknown  boundary  conditions. 

Melt  Phase 

In  the  melt  phase 
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Following  the  same  procedure  outlined  in  Chapter  2 results  in 


pm  _ pml  ^ pm3 

r[0]  ' MO]  + MO] 


+ P 


m4 

[0] 


so  that 


110 


v2pml 

MO] 


= 0 subject  to 


( APml  S 
3P[0] 


m 


S 1 


r=r , 


(E.5.1) 
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Equations  (E.5)  are  solved  using  the  separation  of  variables  technique,  which  results 
in 


+ E OoMcosM 
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(E.6.1) 


(E.6.2) 


, 111 


oo 


[0] 
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(E.6.3) 


and 


OO 


P[0]  = C"’  + £ CTn  1(lanr)coianz) 


(E.7) 


where  an  and  bn  are  given  in  Chapter  2.  The  indicated  forms  of  the  equations  for 
the  pressures  arise  from  the  fact  that  a0=bQ  = 0 and  that 
Iq(0)  = J0(0)  = cos(0)  = cosh(0)  = 0. 

The  constants  are  determined  in  terms  of  the  unknown  boundary  conditions 
in  the  following  manner: 

1)  expand  g 'J  in  terms  of  cos (aj)  and  in  terms  of  J Q(bnr)  so  that 


2)  differentiate  the  pressure  with  respect  to  the  appropriate  variable  for  gx 
and  then  equate  to  the  expansion  of  gx.  Using  as  an  example 
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It  is  apparent  that  g10  = 0 and  since  the  cos^^)  are  independent  of  each  other 
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Performing  similar  calculations  for  g'V  and  results  in 
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Encapsulant  Phase 


In  the  encapsulant 


subject  to 
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Following  the  procedure  described  above  results  in 


(E.9) 


where  cn  and  P0  are  given  in  Chapter  2.  The  pressure  difference  is  given  by  equation 
(E.7)  - equation  (E.8). 

Perturbing  the  right-hand-sides  of  equations  (E.4)  results  in 
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Using  equations  (4.1)  and  (4.2)  in  equations  (E.10),  expanding  in  terms 
of  sin^);  g^  2,  g ™ in  terms  of  cos^^);  g^  4 in  terms  of  J0(^)  and  #3  4 in  terms  of 
P0(c*r)  results  in 
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where 


where 


if  abs  [an-akyai 
if  abs  (an-ak)=ai 


OO 

ftx  m 

S3 k ~ E ^l(0^3kl 


m 

^3 kl  ~ 


+ a 


m 
4 kl 


%4k  ~ E 11 1(0  ^4kl 
1=1 


m - 

y 4kl  - 2bk 


m 

*3 kl 


a 


m 


4kl 


tanh|b  kf3  j cosh  (bkfi} 


OO 


sic  = E ^/(ojy 


e 

1 kl 


(E.11.2) 


(E.11.3) 


(E.11.4) 


where 


117 


e 

y-ikl 


= 2 a, 


hfakTl)  e Kl(«^l) 

T7  r a2w — 7 — r 


ll[ak)  ““  K2^l). 

-E^k-i^-i]^ 


/»=!  ^rl 


+ Xo 


V 


'1/1 


2^  ♦ 2^7-^ 
? n r / i 


1 


+ c 


1 n 


Hanl 


F2 n (rl)  - 2 

+ Jy 

2 n i 


oo 


-E 


(X) 

-E 

/»=1 


al+ak 

+ 

arak 

firt 

(a/+a*)2  + c«  (a/" 

\ 2 2 
fl*J  + C„J 

4i7'  , i 

/ 

al+ak 

^sinhjc^)  2c^ 

\al+ak f 

cl  { ai 

arak 


ai+ak  arak 


where 


4 - E (E'n'5) 


e ~ 2 

ym  - ^k 


118 


oo 

S3 k = £ *l(0?3kl 

where 


(E.11.6) 


+ a 


e 

4 kl 


e 

84kl 


OO 


J2  *1(0?  id 


(E.11.7) 


where 


e o i. 

ym  - 


a 


3kl 


a 


4 kl 


tanh^js)  cosh  (ckp} 


Substituting  equations  (E.ll)  into  equations  (E.7)  and  (E.8)  and  subtracting  results 
in 
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which  is  the  desired  form  of  the  perturbed  pressure  difference. 


(E.12) 


APPENDIX  F 

FOURIER  COEFFICIENTS 


Several  fourier  expansions  are  required  in  this  derivation.  In  general,  a 
function  of  r in  the  melt  will  be  expanded  in  a fourier  bessel  expansion  of  rJ^^y), 
a function  of  r in  the  encapsulant  will  be  expanded  in  a fourier  bessel  expansion  of 
r^i(cnr)'  an^  a z function  in  both  the  melt  and  encapsulant  will  be  expanded  in  a 
fourier  sine  series  of  sin(fly).  Using  K to  illustrate  the  method, 
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where  use  of  the  orthogonality  of  the  functions  has  been  used.  The  integrals  of 
bessel  functions  are  evaluated  through  integration  by  parts.  The  same  procedure  is 
used  for  the  sine  expansions.  All  of  the  integrals  which  arise  from  these  expansions 
are  tabulated  in  standard  references. 

Thus,  the  fourier  coefficients  of  the  F functions  are: 
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